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Abstract: Realistic image synthesis is one of the most relevant subjects in computer 
graphics, and many algorithms have been developed to try to reproduce the visual 
complexity perceived in the real world. However, rendering realistic images in real 
time remains a challenge even with the support of modem graphics hardware. 
Precomputed Radiance Transfer (PRT) is a new graphics technique capable of 
synthesizing highly realistic images in real time. This is achieved by restricting the 
solution of the rendering equation to a subset of the light transport paths that handle 
only energy exchange among diffuse surfaces. Due to the quality of its results, PRT 
has attracted the attention of many computer graphics researchers and practitioners. 
However, understanding and implementing PRT requires familiarity with concepts 
such as projection into basis functions, empirical function integration and light 
transport theory. This tutorial provides a gentle introduction to PRT and its required 
background, enabling the readers to understand and implement the technique. 
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1 Introduction 

Global illumination algorithms are at the heart of photorealistic image synthesis, but 
the high cost usually associated with them has limited their use in real-time applications. 
Because of that, a common practice has been to use algorithms that only implement part of 
the light transport defined by the rendering equation [12], More specifically, a solution for 
the transport among diffuse surfaces in static scenes can be precomputed and then used for 
real-time rendering. This is the case of radiosity [9] methods and, more recently, 
precomputed radiance transfer (PRT) [24], This new class of algorithms is capable of 
producing high-quality renderings, can take advantage of recent developments in 
programmable graphics hardware, and has already been incorporated into Microsoft's 
DirectX SDK [7], However, its use by the research community and practitioners is still timid, 
probably due to the technique's elaborate formulation, which has traditionally involved the 
use of Monte Carlo methods [11] for approximating the solution of some integrals of 
empirical functions, and the use of spherical harmonics as a basis for the space of both 
lighting and transfer functions [6], 

This tutorial presents a step-by-step introduction to precomputed radiance transfer, 
providing a detailed treatment of all aspects necessary for understanding and efficiently 
implementing the technique. In particular, it focus on practical aspects not covered in 
sufficient detail in previous publications. The presentation is illustrated with sample code 
extracted from a real application, which can be used as a reference by researchers, game 
developers and graphics programmers in general. The goal is to provide the readers with a 
solid intuition and understanding of the underlying algorithms, enabling them to incorporate 
the technique in their future applications. 

The tutorial is organized as follows: Section 2 presents a review of the rendering 
equation. Section 3 introduces the notion of precomputed radiance transfer and derives the 
PRT equation from the rendering equation. Sections 4 and 5 provide the background on 
spherical harmonics and Monte Carlo integration, respectively, necessary for understanding 
and implementing PRT. Section 6 describes a PRT implementation, providing code 
fragments for all the core functions. Finally, Section 7 discusses possible extensions to the 
plain PRT rendering, including the incorporation of specular highlights and the use of tone 
mapping. 


2 Illumination Models and the Rendering Equation 

Real world scenes tend to be visually rich (Figure 1), resulting from a complex 
interaction among several factors, such as object geometry, material properties, and the many 
paths followed by the light on its way from its source to the object surfaces. In computer 
graphics, illumination models are abstract representations that attempt to express the 
interplay between light and object surfaces, and can be classified as local and global. Local 
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illumination models compute the shading at a point on a surface considering only its local 
properties. As a result, such a simplification does not handle blocking light {i.e., occlusions), 
which is responsible for producing shadows (an important element for achieving realism), 
and also ignores indirect lighting {i.e., light reflected from other objects that ends up arriving 
at the given point). Global illumination models, in contrast, take into account the entire scene 
when computing the illumination at a surface point. While they can produce more 
sophisticated renderings, this comes at higher computational costs, which has traditionally 
prevented their use in real-time applications. 



Figure 1: A real world scene exhibiting a variety of geometric shapes, 
different materials and light paths. Its appearance results from the 
interplay of all such elements. 


For the specific case of static scenes containing only diffuse surfaces, real-time 
renderings can be achieved using pre-computed solutions for the light transport. This tutorial 
describes one such technique called precomputed radiance transfer. The reader should notice, 
however, that with the advent of modem programmable graphics hardware capable of 
performing several computations in parallel, the gap separating global illumination 
algorithms and real-time performance gets thinner every day. This opens up exciting new 
possibilities and the prospect for rendering more complex light phenomena in real time on 
commodity hardware. 

2.1 The Rendering Equation 

All rendering algorithms can be seen as solutions to particular formulations of a 
general expression known as the rendering equation [12] (Equation 1). Its interpretation can 
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be stated as: the light leaving any point x in a given direction is computed as the amount 
of light that x emits in the direction & 0 , plus the reflected/scattered light from x in the 
outgoing direction to 0 after it has reached x from all incoming directions. 

L(x, coj = L e (x ,co 0 ) + L(x' .Wi) p(x, d>,, coj V{x,x') G{x,x') d cf>, (1) 

where 

Q represents the domain of all possible directions 

L(x', toj is the amount of light arriving at x from another point x ' along <o i 

p (x, fo ; , coj is a function that tells how much of the incoming light 

arriving at x along the direction co,. is reflected along the outgoing 
direction co a 

V(x,x') is a binary visibility function involving x and x' 

G(x, x') is the geometric relationship between x and x' 

In Equation (1) one can identify the factors responsible for the visual complexity 
observed in real scenes: local geometric aspects are implicitly defined by x; material 
properties are expressed as p; the incoming light is represented by L. Occlusions are 
obtained as the result of the visibility function V, so it is related to light transport. Since the 
rendering equation is defined recursively (L appears in both sides of the equation), it is 
useful to think about it as an infinite series, expressed by the Neumann expansion: 

L{x,mJ= L e (x,w 0 )+ L 0 (x,w o ) + L l (x,aj 0 )+L 2 (x,w 0 )+ ... (2) 

The first term of the expansion ( L e ) represents the light emitted by x in the direction 
w 0 . The second term (L 0 ) represents the light reflected/scattered at x in the direction after 
arriving directly from light sources or other emitting surfaces along all possible incoming 
directions. The following term (L } ) represents the light reflected/scattered at x in the 
direction «T 0 after it has bounced once before reaching x. Likewise, L h i>2, represents the 
light reflected/scattered at x in the direction (v a after it has bounced (or been transmitted 
through) i times on scene before reaching x. 

All terms besides L e and L 0 represent indirect lighting (bounces). The first level of 
indirect lighting can be computed using the results obtained with directing lighting. 
Similarly, a second bounce can be performed with the results of the first one, and so on. This 
behavior is expressed by Equation (3), which describes the indirect light contribution exiting 
a point x along the direction <o„ after n bounces. Notice the emittance term has been omitted 
since it has already been accounted for in L 0 . Also, the visibility function must be inverted, 
as one must only consider, for indirect lighting purposes, directions that were blocked when 
direct lighting was computed. 

L n {x,C3 o ) = j Q L (xd3 ( .) p(x, co,.. d>„) (1 - V(x ,x')) G (x, x ') dux (3) 
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Despite its elegance and compactness, directly solving the rendering equation is not 
practical. Its integration domain consists of infinite directions and its recursive nature also 
leads to an infinite number of levels to solve for. Precomputed radiance transfer is an attempt 
to surpass the challenge of real-time global illumination by splitting the problem into two 
parts: one that is solved in a precomputed step and another that is solved in real time. The 
PRT equation can be derived directly from the rendering equation though a series of 
simplifications and approximations. 


3 Precomputed Radiance Transfer 

PRT is based on two main assumptions: (i) all objects in the scene are non-emitters, 
and (ii) the light sources are infinitely distant. The second assumption makes the incoming 
light direction independent of the position of point x. As a result, the term L 0 (direct lighting) 
in Equation (2) is illustrated in Figure 2 and rewritten as: 

L 0 (x, ujj = J fl L € ( dB f ) p(x, to., coj V(x,x') G(x,x') d cu. (4) 

Another common assumption in PRT is to treat all surfaces in the scene as 
Lambertian reflectors. In this case, p becomes just the surface reflectance divided by tt [3], 
and can be taken out of the integral, substantially simplifying the equation. In this tutorial, 
however, we will continue to use the more general form of Equation (4), postponing such a 
simplification until the implementation section. The idea is to give the reader a better 
appreciation for the potential of the technique. The functions P, V and G can be grouped into 
a single function T, known as transfer function. 

L n {x, coj = J fl Lfid,.) T(x,cu o ,w.,x') dw t (5) 

Basically, the lighting function L e describes the arriving light intensity at point x, 
while the transfer function expresses how x responds to the incoming illumination. In PRT, 
these functions are independent from each other and, therefore, can be estimated separately 
and combined later to produce the final result. 




Figure 2: To determine the illumination of a point over a surface, it is necessary 
to integrate, along with its material properties, the lighting function (leftmost), the 
visibility function (center) and the geometric factor (rightmost). 


Even with such simplifications, Equation (5) is not yet practical for real-time 
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implementation. First, one needs a way to estimate the functions L e and T, none of which has 
an analytical solution. And second, we need to efficiently evaluate Equation (5) during 
runtime. In PRT, the lighting and transfer functions are estimated during a preprocessing step 
and their values are used to approximate the integral in real time. These operations will be 
discussed in more detail later. Next, we describe the use of basis functions to project and 
reconstruct arbitrary functions. These operations are essential for understanding how PRT 
can evaluate an estimate of Equation (5) in real time. Figure 3 illustrates how lighting and 
transfer functions can be combined to produce the final shading. 



(a) (b) (c) 



(d) (c) (e) 

Figure 3: The visibility term (a) is combined with the geometric factor (b), which 
gives a transfer factor (c). This transfer factor (c) will be then combined with the 
incoming radiance (d) to determine the final illumination (e). 


3.2 Projection and Reconstruction of Functions 

Given a function space, a basis functions P is an infinite set of functions used to 
project and reconstruct arbitrary functions. The projection of a function / into any function 
Pi gives a measure c, of how similar these two functions are, and is accomplished by 
integrating the product of / and P,over the entire domain of /: 

c i = J /(*) PiW dx ( 6 ) 

The original function / can be recovered as a linear combination of the basis 
functions Pi, each modulated by its corresponding coefficient c,. This process is known as 
reconstruction and its accuracy depends on how many terms are added. The more basis 
functions are used, the better the reconstruction. In fact, when using an infinite number of 
basis functions, which requires an infinite number of coefficients, one is guaranteed to 
recover the original function: 
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f(x) = lim ±c, P,(x) 

By using just a finite number n of terms will produce an approximation / of the 
original function /. 

f(x)~f(x)=£c k P k (x) 

For PRT purposes, the use of orthonormal basis functions, i.e., functions exhibiting 
the following property, is highly desirable: 

JftW £,(*)* =fj wh u en 

[1 when i = j 

Orthonormality guarantees that the integral of the product of two basis functions will 
be either zero, if the functions are different, or one, in case they are the same. Thus, let c k and 
d k be the coefficients associated with the projections of any two functions /( x ) and g{x% 
respectively, into the basis function P k , for all k. The integral of the product f{x)g(x) can 
then be obtained as the sum of the products c k d k . More formally: 

c k = ! f(x) p k (x)dx 

d k =fg(x)p t (x)dx 

J fix) g(x) dx « J Z c k P k {x) Z d k p k (x) dx = 

Z c k d k J P k {x) P k {x) dx = 

Z c k d k (7) 

This provides a powerful mechanism for factoring the evaluation of Equation (5) in 
two steps: (i) the projection of both L e and T into an orthonormal basis, which can be done 
during preprocessing, and (ii) the evaluation of the integral as a dot product (Equation 7), 
which can be efficiently done in real time. This will be explored in detail next. 

3.3 Precomputed Radiance Transfer Equation Derivation 

Let both lighting and transfer functions from Equation (5) be projected into the same 
set basis functions fi k . In this case, the original functions can be approximated as: 

PkWlll where /*=/ L € [w 4 ) P k {w t )d0 t (8) 
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and 

T(x,io 0 ,(b i , x') » X h where t k =f T(x,uj 0 ,co i , x') jS t (u5 ; ) dw i (9) 

Substituting the value of the approximate lighting function (Equation 8) into 
Equation (5): 

L 0 (x, coj = J Q [ X h Pk(&j) T(x,co 0 ,a> jf x') ] dc5 i 

Using the fact that integration is a linear operation - the integral of sums equals the 
sum of integrals - the equation above can be rewritten as: 

L 0 (x, cuj = X h J fl Pk(^i) T{x,w 0 ,cb i ,x')du) i 

For each term k of the sum, the associated lighting coefficient will be multiplied by 
the entire integral. But this integral represents a projection operation. So, each term of the 
sum will produce a new coefficient that is the result of projecting the transfer function into 
the associated basis function, resulting in the following equation: 

Hix,dS 0 ) = il k t k ( 10 ) 

This equation is identical to Equation (7), and it can be seen as a dot product of the 
lighting and the transfer coefficients. One should note that Equation (10) only handles direct 
lighting (i.e., the light leaving x after having arrived there directly from the light source). 
Notice, however, that it handles occlusions with respect to the light source. To be able to 
perform extra bounces, a similar derivation can be done with higher order terms of the 
Neumann Expansion, which leads to the final PRT equation, shown below: 

L(x, «3J = ^ l k {t° k + t\ + t\ + ...) (11) 

It can be seen that one set of transfer coefficients should be obtained for each 
outgoing direction, which is explicit in the transfer function T. This is the general case, 
which can be simplified to a single set of transfer coefficients (i.e., independent of the 
outgoing direction) by assuming that all surfaces are ideal diffuse reflectors (Lambertian 
surfaces). Such a simplification will be discussed later in Section 6. 

At this point, one is still left with two important questions: which basis function 
should be used and how to solve the integrals required to perform the projection? These 
questions will be addressed in the next sections. 
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4 Spherical Harmonics 

Spherical harmonics (SH) are orthonormal functions defined over the unit sphere, 
where the 2D domain can be seen as the set of all possible directions. In their general 
formulation, SH functions are defined over complex numbers, but for PRT purposes, we are 
interested in approximating real functions. Therefore, this tutorial will focus on real spherical 
harmonics only. 


4.1 The Legendre Polynomials 


Legendre Polynomials are the core of the spherical harmonic functions. The 
Associated Legendre polynomials return real numbers and are recursively defined as follows: 


{/ elN } 

{ nzeIN | m<l } 

{xeIR | -1<x<1 } 




x{2m+l )P m m 

(-l) ra (2m —1]/1 (l—x 2 ) m/2 
x(2l-l)PJL I -(l + m-l)PJ!% 


(l-m) 


when l=m + \ 
when l = m 

otherwise 


They are characterized by two parameters, /, which is usually called band, and m, that varies 
according to /, and they are defined over real numbers in the interval between [—/,/]. The 
unusual operator / / is called double factorial, and is defined as follows: 


' (»(*- 2)1 


when n< 1 
when n> 1 


4.2 Spherical Harmonic Functions 

Spherical harmonic functions are parameterized using 9 and 4>, as shown below, with 
scaling factors K"‘ for normalizing the functions. 

{/ elN ) 

(meZ | —l<m<l } 

[V2 K™ cos{m4>) PT(cos( 9)) when m >0 
T/ (9, (/>) = W2 K” sin(-m4>) P^icosfi)) when m< 0 
P°,(cos{9)) when m = 0 

where 
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l (2/+l)(/-|iw|). 

4tt(/+M) 


A spherical harmonics approximation of a function using l bands requires / 2 
coefficients. An alternative representation using a single index can be obtained using the 
following relationship: 


{y”{6,4>)=y i {6,ct>) \ i=l(l+l)+m] 


That is all one needs to know about spherical harmonics to use them in PRT. The next 
section discusses Monte Carlo integration, a method for numerically evaluating the integral 
of any function, be it analytical or empirical. 


5 Monte Carlo Integration 

The expected value E of a function / with a random variable x associated to it is an 
average value that will tend to return most often if evaluated to a large number of samples. A 
random variable is a value that lies within a specific domain and has a probability p 
associated to it. Thus, E can be obtained by computing the following integral over the entire 
domain of the random variable: 

E[f{x)\=jf(x)p{x)dx (12) 

Another way to get the expected value of a function is to take the mean of an infinite 
number of random samples within its domain: 

E[f(x)] = limiX/(x ; ) (13) 

Using a finite number of random samples will give only an approximation of the 
expected value: 

(14) 

Combining Equations (12) and (14), one gets a numerical solution for estimating the 
integral of an arbitrary function, which is known as Monte Carlo integration (Equation 15). 

j /WA = j/MzW*.I|;fg| (15| 

To solve an integral using the Monte Carlo, it is then necessary to take lots of samples 
of the function and each one must have an associated probability. Equation (15) can be 
rewritten in terms of a weight function w as: 
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J/(*)<fc~±Z/(*,)*(*,) 


(16) 


5.2 Using Monte Carlo to Solve the Projection into Spherical Harmonic basis 

Replacing the general basis functions P used in Equation (6) by the spherical 
harmonic basis functions y and parameterizing the function / in terms of the space of 
directions, Equation (6) can be rewritten as: 

c?= f f(w) y”(w) dw (17) 

and it is useful to think of both spherical harmonic functions and coefficients in terms of a 
single index. In this case 

c k = $ f(&) yk(“>) dw (18) 

Applying Monte Carlo integration, the coefficients obtained with the projection 
operation can be approximated using n samples as shown in Equation (16). Note that the 
more samples are used, the better the approximation. It should be emphasized that each 
sample /(co,-) is associated with a direction cSj. Thus, both the lighting and transfer 
functions need to be sampled along the set of directions m,, 1 < j<n. 

c k M ^ Z/(«5y) yniwjwi&j) (19) 


Although one knows how to evaluate both the function / and the SH basis functions 
for an arbitrary direction, it is still necessary to associate a probability to the occurrence of 
this direction. Since these functions are defined over the space of all possible directions 
(geometrically, this space can be seen as the surface of a unit sphere), the probability of any 
sample (direction) occurring in such a surface will be one over the area of this unit sphere 
(i.e., 1 /4tt). So, the probability associated to each directions is: 


P(&j) = — 


Substituting the weight function in Equation (19), one gets Equation (20), which 
describes the projection of an arbitrary function / defined over the space of all possible 
directions into the spherical harmonic basis. 

Z f(&j) ?*(«*;) ( 2 °) 

Since our Monte Carlo formulation assumes that each sampling direction is equally 
probable 1 , the actual samples of the function / should be evenly distributed over the unit 


1 An alternative to this approach is the use of importance sampling. 
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sphere. Moreover, the same set of sampling directions m, should be used for both the 
lighting and the transfer functions. In order to guarantee a uniform sampling of the unit 
sphere, the following algorithm known as stratified sampling is often used: 

(i) Evenly distribute the number n of samples over the unit square. For this, subdivide 
the unit square into Vn X 4n cells, and randomly select a sample inside each cell; 

(ii) Map the coordinates of the samples in the unit square to coordinates on the unit 
sphere using Equation (21). 

(x, >•) -* ( 2 acosU\-x), 2ny) -> {0, f) (21) 

The idea of stratified sampling is illustrated in Figure 4. On the left, one sees the 
distribution of the generated samples on the unit square. Each sample on the square has been 
mapped to one sample on the sphere on the right, according to Equation (21). 



Figure 4: Uniform sampling of the sphere using stratified 
sampling. A total of 10000 random samples generated over the 
unit square (left). These samples are mapped to the surface of 
the sphere on the right using Equation (21). 


6 PRT Implementation 

At this point, the reader has all the background necessary to implement PRT, which is 
the subject of this section. As mentioned before, the PRT algorithm consists of two steps: (i) 
the projection of the lighting and transfer functions into some orthonormal basis functions, 
which happens as part of preprocessing, and (ii) the evaluation of Equation (11), which is 
done in real time. 

For this tutorial, lighting functions will be represented as environment maps stored as 
spherical maps (light probes), but it is also possible to use cube maps. We also briefly 
discuss the use of analytical representations for lighting functions. For the case of transfer 
functions, two types will be considered: (i) unshadowed transfer functions, which ignore the 
visibility function (V(x,x ') in Equation 4), and (ii) shadowed transfer functions, which take 
the occlusion term into account. We also briefly discuss what it takes to add indirect lighting 
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(also called interreflected transfer function) to the shadowed case. The transfer function will 
be evaluated for each vertex present in the scene. 

The code fragments shown in this section were written in C/C++. There is also 
sample code for GPU execution written in Cg [15], with its corresponding CPU version. 

6.2 Precomputation Step 

This is the core of a PRT implementation, concentrating most of the coding effort. 
The following sections discuss each of the precomputation steps. 

6.2.1 Sampling Directions 


The first step is to generate the m sampling directions using the stratified sampling 
strategy described in Section 5.2 and illustrated in Figure 4. Before generating these 
directions, it is necessary to define the sample data structure: 


struct 

: vector3 


floa 

t x; 


floa' 

t y; 


floa 

}; 

t z; 


struct 

: Spherical 


^ f'1 oa 

t theta; 


floa 

} 

t phi; 


struct 

: Sample 


Sphe 

rical spherics 

il_coord; 

vecti 

or3 cartesian. 

.coord; 

floa 

}; 

t* sh_functior 

is; 


Each sample stores its spherical coordinates ( 6,4 >) and its Cartesian coordinates 
(x,y,z). The spherical coordinates will be used to evaluate the spherical harmonic 
functions, while the Cartesian ones will be used to evaluate the lighting function (compute 
dot products). The samples also store the result of the evaluation of the specified bands of the 
spherical harmonic functions (for the associated spherical coordinates). The details behind 
this extra storage will be explained later, 
struct Sampler 

Sample* samples; 
int number_of_sampfes; 


The sampler data structure (above) stores the samples, providing access and the 
ability to iterate over the samples (created using stratified sampling). The conversion 
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between spherical coordinates obtained with Equation (21) to Cartesian coordinates is given 
by: 


( 9, cf>) —> (sin0cosd>, sin0sind>, cos0) —> (x,y,z) 

The following code fragment performs the sampling and initialization of the Sampl er 
data structure with nxn samples. The sh_functions variable is initialized with null and, in 
a later step, memory will be allocated to store the values for the evaluated bands. 


void GenerateSamples(Sampler* sampler, int N) 

Sample* samples = new Sample [N*N] ; 
sampler->samples = samples; 
sampler->number_of_samples = N*N; 
for (int i =0; i < N; i++) 

for (int j = 0; j < N; j++) 

float a = ((float) i) + Random()) / (float) N; 

float b = ((float) j) + Random()) / (float) N; 

float theta = 2*acos(sqrt(l-a)); 

float phi = 2*Pi*b; 

float x = sin(theta)*cos(phi); 

float y = sin(theta)*sin(phi); 

float z = cos(theta); 

int k = i*N + j ; 

sampler->samples[k].spherical_coord.theta = theta; 
sampl er->sampl es [k] . spherical_coord. phi = phi; 
sampler->samples[k].cartesian_coord.x = x 
sampler->samples[k].cartesian_coord.y = y 
sampler->samples[k].cartesian_coord.z = z 
sampler->samples[k].sh_functions = NULL; 


} 

}; 

float Random() 

{ float random = (float) (rand() % 1000) / 1000.Of; 
return(random); 

} 


As it is a C based implementation, the rand() function works properly, a random seed 
must be initialized. So, at the early stages of the main() function, it can be initialized with 
the call showed below: 
s rand(time(NULL)); 


6.2.2 Spherical Harmonics Evaluation 

In order to evaluate real spherical harmonic functions, one needs to be able to 
evaluate the associated Legendre polynomials, which can be done with the following code 
fragment. 
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float Legend re (int 1, int m, float x) 

float result; 
if (1 == m+1) 

result = x*(2*m + l)*Legendre(m, m); 
else if (1 == m) 

result = pow(-l, m)*DoubleFactorial (2*m-l)*pow((l-> 
else 


} 

float DoubleFactorial (int r 

if (n <= 1) 
return(l) ; 
else 

return(n * DoubleFactori 

} 


Now, the spherical harmonic functions can be evaluated 


float SphericalHarmom 

float result; 
if (m > 0) 
result = sqrt(2) * 
else if (m < 0) 
result = sqrt(2) * 
else 


c(int 1, int m, float theta, float phi) 

K(1, m) * cos(m*phi) * Legendrefl, m, cos(theta)); 
K(1, m) * sin(-m*phi) * Legendre(l, -m, cos(theta)); 
Legendre(l, 0, cos(theta)); 


where the normalization factors can be computed as: 


float K(int 1, int m) 

{ 

float num = (2*1+1) * factorial (l-abs(m)); 
float denom = 4 *pi * factorial (l+abs(m)); 
float result = sqrt(num/denom); 
return(result); 


The presented implementations for both the associated Legendre polynomials and for 
the double factorial are recursive. They were used here for simplicity, as the mappings 
between their mathematical definitions and implementations are straightforward. Also, these 
functions are only used in the preprocessing stage. More efficient ways to implement these 
functions iteratively can be found in [20], 


6.2.3 Precomputing the Spherical Harmonic Bands 

It is useful to compute and store the SH bands at the samples to avoid evaluating them 
several times in further steps. The following code fragment does this. 
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void PrecomputeSHFunctions(Sampler* sampler, int bands) 
for (int i =0; i < sampIer->number_of_samp1es; i++) 

float* sh_functions = new float [bands*bands]; 
sampler->samples[i].sh_functions = sh_functions; 
float theta = sampler->samples[i].spherical_coord.theta; 
float phi = sampler->samples[i] .spherical_coord.phi; 
for (int 1 = 0; 1 < bands; 1++) 
for (int m = -1; m <= 1; m++) 

{ int ’ - 1* 1+1 

sh_functions[j] = SphericalHarmonic(l, m, theta, phi); 


This function takes two parameters: the sampler, which was created previously, and 
the number of spherical harmonic bands to be used (the more bands, the better the 
approximation). Recall that an approximation using / bands uses / 2 spherical harmonic 
functions. The code iterates over the bands and evaluates these functions. This process must 
be repeated for all samples. 


6.2.4 Lighting Function Representation 

There are essentially two ways to represent a lighting function: empirically and 
analytically. In this tutorial, light functions are represented empirically using light probes 
(Figure 5). A light probe is an omnidirectional image that records the incident illumination at 
a particular point in space and can be captured from real world scenes. There are several 
freely available light probe images on the Internet [14]. 



Since a light probe is only a 2D image, one needs to convert spherical or Cartesian 


146 


RITA • Volume XIII • Numero 2 • 2006 




A Gentle Introduction to Precomputed Radiance Transfer 


coordinates into image coordinates. This can be done using the following code fragment 
available from Paul Debevec's IBL Tutorial [5], Basically, the function converts the 3D 
Cartesian coordinates (direction) into 2D texture coordinates (tex_coord) and, finally, 
converts them to the pixel offset in the image (pi xel_coord), based on image dimensions. 
The output parameter color will have the RGB color associated to that direction, after its 
execution. 


void LightProbeAccess(Color* color, image* image, vector3 direction) 


{ 


d = sqrtCdirection.x* 
float r = (d == 0) ? O.Of : 
float tex_coord [2]; 
tex_coord[0] = 0.5f + directi on.x * r; 
tex_coord[l] = 0.5f + directi on.y * r; 
int pixel_coord [2]; 

pixel_coord[0] = tex_coord[0] * image.width; 
pixel_coord[l] = tex_coord[l] * image.height; 
int pixel_index = pixel_coord[1] *image.width h 
color->r = image.pixel [pixel_index] [0]; 
color->g = image.pixel [pixel_index] [1] ; 
color->D = image.pixel [pixel_index] [2] ; 


pixel_coord[0]; 


The lighting function can also be represented as a cube map, but it is a lot easier to 
fetch the desired value from a spherical map, such as a light probe. Another way to evaluate 
a lighting function is to use an analytical representation. An example of such a function is 
shown below, which corresponds to two monochromatic light sources, at 90 degrees from 
each other [10], Figure 6 provides an illustration for this lighting function. 

light{e, 0) = max( 0, 5 cos(0)-4) + mox(O, -4sin(0-Tr) cos(</>-2.5)-3) 



Figure 6: A 3D plotting of an analytical spherical lighting 
function that defines two monochromatic light sources, at 90 
degrees from each other. 
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6.2.5 Projection into Spherical Harmonic basis using Monte Carlo Integration 

Everything that is necessary to implement the projection operator is now available: 
the samples, the function, and the SH basis functions. It is time to take a look at Monte Carlo 
integration in action, and it is really simple to implement it. The weight function is just a 
constant and the number of samples is known. The following code fragment implements the 
projection operation using Monte Carlo integration: 


nd ProjectLightFunctioi 


: bands*bands; i++) 


coeff s [i ]. i 
coeffs[i] .< 
coef f s [i ]. f 


= O.Of 
= O.Of 
= O.Of 


for (int i =0; a < 

{ 


vector3& directioi 
for (int j = 0; j 


sampl er->nutnber_of_sampl es; i ++) 

= sampler->sampl es[i].cartesian_coord; 
< bands*bands; i++) 


} 


} 


Color color; 

LightProbeAccess(&color, fight, &direction); 
float sh_function = sampler->samples[i] .sh_functi< 
coeffs[j].r += (color.r * sh_function); 
coeffs[]].g += (color.g * sh_function); 
coeffsEjl.b += (color.b * sh_function); 


;[j]; 


float weight = 4.0f*Pi; 

float scale = weight / sampler->number_of_sampfes; 
for (int i =0; i < bands*bands; i++) 

coeffs[i].r *= scale; 
coeffs[i].g *= scale; 
coeffs[i].b *= scale; 


The code initializes all coefficients with zero. The sum operates over all samples. For 
each sample, the directional Cartesian coordinates are used to access the light probe and the 
spherical harmonic precomputed values associated with each band are retrieved. They are 
combined and accumulated. At the end, all values are scaled using the weight function. 

One important thing to note is that three sets of coefficients were generated 
independently, one for each color channel. For monochromatic lighting functions, one vector 
of coefficients is enough. The same principle applies when projecting transfer functions, 
which will be discussed next. 


6.2.6 Transfer Functions 

Recall that the transfer function, as defined within Equation (5), is the product of 
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three other functions: the surface scattering function p, a visibility function V, and a 
geometric term G that expresses the relationship between the emitter and the receiving point. 
These intuitively simple functions are surprisingly the most difficult to implement in the 
whole PRT technique. 

The function p can be as complex as we want, handling phenomena such as caustics, 
refraction, subsurface scattering, and so on. For this tutorial, a simple scattering function will 
be used, assuming that the surfaces are ideal diffusers (i.e., Lambertian surfaces), meaning 
that they reflect the incident light with equal intensity in all directions. Using such a 
simplification, the dependence between the transfer function and an outgoing direction 
disappears. Also, the scattering function will be a constant value given by [3]: 

P (x , CO:, «3 ) = — 

t r 

The P d term in the equation above is a reflectance coefficient associated to the point x 
where the transfer function will be evaluated, just like the diffuse coefficient used in the 
Phong Illumination Model [19]. 

For the geometric term G, one can use the one known as Lambert's cosine law 
(Equation 22). Here, x ' is the emitter, co, is the normalized incident light direction, and n x is 
the normalized surface normal at x. For some incident directions, the dot product in Equation 
(22) can be negative, representing light arriving from behind the surface. Equation (22) 
clamps such values to zero. 

G(x,x') = max{0,n x -m i ) (22) 

By blocking some light paths, the visibility function is responsible for a great deal of 
realism. Despite its simple definition - returns 0, if the light is blocked; 1, otherwise - it is an 
intricate function to implement. In fact, the visibility function is implemented as a ray¬ 
casting procedure. Before going any further into the details about visibility function 
implementation, it is worthy to take a look at some kinds of transfer function variations. 

6.2.6.1 Unshadowed Transfer Function 

By simply ignoring the visibility function, one can produce renderings similar to that 
obtained with the use of local illumination models, such as Phong and Cook-Torrence 
models [19] [4], but with one advantage: the use of any number of light sources without any 
additional rendering cost. Figure 7 exemplifies how unshadowed transfer works: 

Before implementing the transfer function, it is necessary to have a data structure to 
handle the scene information (e.g., objects, vertices, triangles, materials). The following data 
structure will be used to describe the scene: 
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Figure 7; For the unshadowed transfer function, occlusions are 
ignored. As a result, d5j will be considered when illumination is 
evaluated for x. 


truct Scene 

vectorB* vertices; 
vectorB* normals; 
int* material; 

Triangle* triangles; 
lor* albedo; 

: number_of_vertices; 


Consider a scene composed of triangles, each with three indices to vertices and the 
same indices also used for the normals. Each vertex has an associated material, which is an 
index to a color (albedo). 


for (int i =0; i < 
for (int j = 0; j 


scene->number_of_vertic< 
< bands*bands; j++) 


M\ 


^ coeffs[i] [j] .6 = O.Of 

for (int i =0; i < scene->number_of_vertices; i++) 
for (int j = 0; j < sampler->number_of_samples; j++) 

Sample* sample = sampler->samples[j] ; 

float cosine_term = dot(&scene->normals[i], &sample.cartesiai 
for (int k = 0; k < bands*bands; k++) 

float sh_function = sample.sh_functions[k]; 
int materia_idx = scene->material [i]j; 

Color* albedo = scene->albedo[materia_idx]; 
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coeffs[i] [k]. 
coeffs[i] [k] .1 
coeffs[i] [k]. i 


+= (albedo.r * sh_function * cosine_term); 
+= (albedo.g * sh_function * cosine_term); 
+= (albedo.b * sh_function * cosine_term); 


} 

float weight = 4.0f*Pi; 

float scale = weight / sampler->number_of_samples; 
for (int i = 0; i < scene->number_of_vertices; i++) 
for (int j = 0; j < bands*bands; j++) 

coeffs[i] [j]. r *= scale; 
coeffs[i] []] .g *= scale; 
coeffs[i] Lj] -b *= scale; 


The projection of unshadowed transfer functions into SH basis is pretty similar to 
light projection. The main difference is that each vertex will have a transfer vector (thus the 
double indirection Color** is necessary for the coeffs variable). Inside the sum, the 
geometric factor must be evaluated, which is the dot product of the vertex normal and the 
sample direction. Since each vertex may have a different material, this must be retrieved for 
each vertex. The associated spherical harmonic function is retrieved just as it was for light 
projection. At the end, the resulting value is scaled as before. 



Figure 8: A Teacup rendered using an unshadowed transfer function under three different lighting 
functions. 


The previous code fragment implements the projection of the transfer function into 
SH basis. The transfer function is evaluated using the scattering function (constant) and 
Lambert's cosine law (Equation 22) for the geometric term. Obviously, the visibility function 
is omitted. It is important to notice that the transfer function will be evaluated and projected 
for every vertex in the model, and each vertex will have as many transfer coefficients as the 
number of basis functions used. Recall that this number is the square of the number of SH 
bands used, just like when projecting the lighting function. 

Figure 8 illustrates the use of the unshadowed transfer function for rendering a teacup 
using three different lighting functions. Note the influence of the lighting function in the 
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perceived color of the object. 

6.2.6.2 Shadowed Transfer Function 


In order to produce a shadowed transfer function, the visibility term must be taken 
into account, which can be implemented using ray-casting. The code fragment below 
implements the Moller et al. algorithm for ray-triangle intersection [16]: 


bool 

RaylntersectsTriangle(Vector3* p, Ve 

:ctor3* d, 


vector3* vO, vector3* vl, vector3* v2) 

flc 

>at el [3] = { vl->x - v0->x, vl->y - 

v0->y, vl->z - v0->z }; 

f 1 c 

>at e2 [3] = { v2->x - v0->x, v2->y - 

v0->y, v2->z - v0->z }; 

flc 

•at h [31; 



>ss(h, d, e2) ; 


flc 

>at a = dot(el, h); 


if 

(a > -O.OOOOlf && a < O.OOOOlf) 



eturn(fal se); 


flc 

•at f = l.Of / a; 


flc 

iat s [3] = { p->x - v0->x, p->y - vO- 

->y, p->z - v0->z }; 

flc 

•at u = f * dot(s, h) ; 


if 

(u < O.Of I I u > l.Of) 



eturn(false); 


flc 

iat q [3]; 


crc 

>ss(q, s, el); 


flc 

•at v = f * dot(d, q); 


if 

(v < O.Of || u + v > l.Of) 


return(false); 


flc 

•at t = dot(e2, q)*f; 


if 

(t < O.Of) 


return(false); 


> ret 

:urn(true); 



This code is called inside the visibility function. Since this is an unoptimized code, it 
will check all other triangles of the mesh against the ray. If no triangle is intercepted, the 
direction is free from obstacles, and the visibility function returns 1; otherwise, it returns 
zero. The code fragment below implements the visibility function, 
bool visibility(Scene* scene, int vertexidx, vector3* direction) 
bool visible (true); 

vector3& p = scene->vertices[vertexidx]; 

for (int i = 0; i < scene->number_of_triangles; i++) 

Triangle& t = scene->triangles[i]; 

if ((vertexidx != t.a) && (vertexidx != t.b) && (vertexidx != t.c)) 

vector3& vO = scene->vertices[t.a]; 
vector3& vl = scene->vertices[t.b]; 
vector3& v2 = scene->vertices[t.c]; 

visible = !RaylntersectsTriangle(&p, direction, &v0, &vl, &v2) ; 
if (!visible) 
break; 

} 

} 

return(visible); 
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Figure 9: For the shadowed transfer function, occlusions are 
considered, producing shadows. The direction dl, has been 
occluded by other elements of the scene and will not be 
considered when illuminating x. 


For the case of the shadowed transfer function (Figure 9), the projection code is 
identical to that presented for unshadowed transfer, except for the addition of the binary 
visibility function. 


{ 


Pro]ectshadowed(Color* 
Scene* 


coeffs, 
scene, i 


, Sampler* sampler 
int bands) 


for (int i =0; i < scene->number_of_vertices; i++) 
for (int j = 0; j < sampler->number_of_samples; j++) 

Sample& sample = sampler->samples[j] ; 
if (visibility(scene, i, &sample.cartesian_coord)) 

float cosine_term = dot(&scene->normals[i], &sample.cartesian_coord); 
for (int k = 0; k < bands*bands; k++) 

{ 

float sh_function = sample.sh_functions[k]; 
int materia_idx = scene->material [i]; 

Colors albedo = scene->albedo[materia_idx]; 
coeffs[i][k].r += (albedo.r * sh_function * cosine_term); 

coeffs [i ] [k] .q += (albedo, q * sh_function * cosine_term); 

coeffs [i ] [k] .6 += (albedo. 6 * sh_function * cosine_term); 


} 


1 


Figure 10 shows the result of using the shadowed transfer function for rendering the 
teacup from the same viewpoint shown in Figure 8 and under the same lighting functions. 
Note the considerable increase in realism due to the presence of soft shadows, when 
compared to the images shown in Figure 8. 
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Figure 10: The teacup rendered using a shadowed transfer function under three different lighting 
functions. Notice the soft shadows that greatly increases the realism of the scene. 


One should notice, however, that the plain shadowed transfer function only computes 
one level of transport (i.e., the direct lighting transport level). For additional bounces, some 
considerations should be done, and will be discussed in the next section. 

6.2.6.3 Interreflected Transfer Function 

Given the shadowed transfer function, it is possible to consider some refinements to 
perform indirect lighting and make the renderings even more realistic. However, one should 
note that any additional level of indirect lighting implies more transfer coefficients that need 
to be precomputed, stored and then handled during runtime. However, as it was shown 
previously, these transfer coefficients can be used to compute the next bounce (Equation 3) 
and, when all bounces were performed, they can be merged, resulting in a final transfer 
vector (Equation 2), which handles all bounces plus the one for direct transfer, obtained from 
the shadowed transfer vector. 

To implement interreflected transfer, first the transfer coefficients for shadowed 
transfer should be computed. By doing so, each vertex of the scene will have a transfer 
vector associated to it. The next step is to perform the additional bounces, obtaining the 
terms step by step and adding them at the end. 

Even if it is easy to understand the basics, implementing such a transfer function may 
be challenging (and complicated to debug). But it is just a matter of brute force. For more 
details on how to implement the interreflected transfer function, please see [10]. 

6.3 Real Time Step 

Once computed, the light and transfer coefficients are ready to be used for real-time 
rendering. This is the easiest part of the entire technique, which consists of computing a per- 
vertex dot product — one that can have lots of coefficients, but still a dot product. 
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6.3.1 Evaluating the PRT Equation in the CPU 

To determine the final shading of a point (vertex), it is necessary to evaluate the PRT 
equation (Equation 8 or 9, depending on the existence of indirect lighting) for it by 
computing the dot product between its transfer vector and the lighting vector. The following 
code fragment computes the dot product for each color channel and renders the results using 
OpenGL (with Gouraud shading). 


oid Render(Color* light, Color** coeffs, Scene* scene, 
glBegin(GL_ triangles); 

for (int i = 0; i < scene->number_of_triangles; i++) 

Triangle& t = &scene->triangles[i]; 
vector3& vO = scene->vertices[t.a]; 
vector3& vl = scene->vertices[t.b]; 
vector3& v2 = scene->vertices[t.c]; 


Color cO = { O.Of, 
Color cl = { O.Of, 
Color c2 = { O.Of, 
for (int k = 0; k 

cO. r += (light 
cO.g += (light 
cO.b += (light 
cl. r += (light 
cl.g += (light 
cl.b += (light 
c2.r += (light 
c2.g += (light 
c2.b += (light 


O.Of, O.Of }; 

O.Of, O.Of }; 

O.Of, O.Of }; 

: bands*bands; k++) 

r * coeffs 
g * coeffs 
b * coeffs 
r * coeffs 
g * coeffs 
b * coeffs 
r * coeffs 
g * coeffs 
b * coeffs 


■ r) 

■9? 


■ o ) 

■ r) 

■g? 


glColor3f(c0.r, cO.g, cO.b); 
glvertex3f(vO.x, vO.y, vO.z); 

glColor3f(cl.r, cl.g, cl.b); 
glvertex3f(vl.x, vl.y, vl.z); 


} 


glColor3f(c2.r, c2.g, c2.b); 
glvertex3f(v2.x, v2.y, v2.z); 


} 


gl End(); 


bands) 


6.3.2 Evaluating the PRT Equation in the GPU 

Evaluating the PRT equation in the GPU is pretty straightforward. All it takes is to 
declare one uniform parameter ( i.e ., a parameter that will remain constant over the 
processing of all vertices), which contains the lighting function coefficients. It is still 
necessary to declare the transfer vector, which varies for each vertex. Given that, the vertex 
program only needs to compute the dot product and forward the result to the rasterizer. 
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The following code implements the vertex shader for PRT: 

struct app2vertex 


float4 f4Position 

float4 f4Color 

float3 vf3Transfer [N*N]; 

}; 

: POSITION; 

: COLOR; 

struct vertex2fragment 


float4 f4ProjPos 
float4 f4Color 
}; 

: POSITION; 

: COLOR; 

vertex2fragment VertexShader 


app2vertex IN, 

uniform float3 vf3Light [N*N], 
uniform float4x4 mxModelviewProj 

{ 

vertex2fragment OUT; 


OUT.f4ProjPos = muI(mxModel 

ViewProj, lN.f4Position); 

OUT.f4Color = float4(0.0f, 
for (int i =0; i < n*n; i4 

O.Of, O.Of, 1.Of); 

+) 

{ 

OUT.f4Col or.r += (lN.vf3T 
OUT.f4Color.g += (lN.vf3T 
OUT.f4Color.b += (in. vf3 t 

} 

ransfer[i] . r * vf3Light[i]. r) ; 
ransfer[i].g * vf3Light[i] .g) ; 
ransfer[i] .6 * vf3Light[i] .6); 

return(OUT); 

} 



The following code implements the pixel shader for PRT: 


struct fragment2screen 


float4 f4Color 
}; 

: COLOR; 

vertex2fragment Pixel Shader 


vertex2fragment IN 


{ 

fragment2screen OUT; 

OUT.f4Color = lN.f4CoIor; 
return(OUT); 

} 



7 Extensions to the Plain PRT Rendering 

PRT provides a global illumination solution for scenes containing only diffuse 
surfaces. To increase the realism of the generated images, it is possible to add specular 
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effects to the final rendering. Sloan et al. [24] described a method to render glossy reflections 
using PRT. It is also possible to add high frequency specular reflections using reflection 
mapping [2] [8], In this case, the results will not be physically correct, but the technique is 
very easy to implement, very fast to evaluate in real time, and produces more sophisticated 
images. Figure 11 shows the result of integrating diffuse PRT with reflection mapping. On 
the left, pure diffuse and specular renderings are shown. On the right, both techniques have 
been combined, producing a more realistic image. A code fragment showing how to combine 
the diffuse and specular contributions is shown below. 

The reflection mapping implementation should use an environment mapping that is 
consistent with lighting function (light probe). To implement reflection mapping on the CPU, 
one can use the EXT_texture_cube_map OpenGL extension, for which a tutorial can be 
found at [18]. For a GPU implementation of reflection mapping, refer to the NVIDIA SDK 

[17]- 


When using high dynamic range images for the environment texture, a blooming 
effect [17] and tone mapping [22] [27] can also be used to enhance the the image quality 
even further. These techniques were used to synthesize Figure 11 (right). A comprehensive 
discussion on how to combine diffuse PRT, reflection mapping and high dynamic range 
images in real time can be found in [26], 


diffuse 

diffuse + specular 


material specularity: 
teacup: 0.3 
plate: 0.1 

specular 

*' . \ 

* Wm 

^ _^ N 

W - - ^ , 



Figure 11: A teacup rendered using both diffuse PRT and reflection mapping. An high dynamic 
range image was used for lighting function and environment texture, so tone mapping was 
employed to properly exhibits the generated image. 
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Once both diffuse PRT result and reflection mapping color have been computed, one 
can combine them using the following formula, where specular factor is equivalent to the 
specular coefficient used in the Phong model [19]. Remember that the diffuse coefficient was 
already encoded in the PRT diffuse color, as explained on section 6.2.6: 

FinalColor =PRTDiffuseColor + SpecularFactor* Specular Color 

8 Conclusion 

Precomputed radiance transfer is capable of rendering highly realistic images of static 
diffuse environments in real time. It assumes an infinitely distant light source, and consists of 
factoring the rendering equation into incident light and light transport terms. By projecting 
these terms into a family of orthonormal basis functions, the resulting simplification of the 
rendering equation can be estimated using dot products. PRT has a big advantage over the 
commonly used local illumination models since it is able to render scenes using any number 
of light sources without compromising performance. It is also capable of producing soft 
shadows that highly increases the realism of the scenes. By taking extra light bounces into 
account, one can simulate diffuse interreflection among the objects in the scene, increasing 
the realism even further. 

PRT can produce glossy reflections [24], and specular ones can be easily combined 
with PRT renderings by using reflection mapping. In this case, the resulting images tend to 
look quite realistic even though the approach is not physically correct. 

Currently, PRT cannot be used with dynamic scenes. Although, spherical harmonic 
functions allow the lighting function or the scene to be rotated, preserving the precomputed 
light transport, it does not allow object translation, since this would change the visibility 
function. Some research has been done to try to remove this limitation [1] [13] [23] [25] 
[28], but no definitive solution has been found yet. 

The quality of the renderings obtained with PRT depends on the basis functions 
chosen to perform the projections. Using a small number of coefficients, spherical harmonic 
functions are acceptable just to approximate low frequency characteristics of the lighting 
function. One way to support high-frequency approximation of the lighting function with less 
coefficients is to use Haar Wavelets [21], 

This tutorial focused on a gentle introduction to PRT, using diffuse reflection and 
direct lighting only. More sophisticated transfer functions are possible supporting, for 
instance, interreflections, caustics, refractions, and subsurface scattering. PRT has proved to 
be a very interesting and promising rendering technique. 
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